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Preface 


The  Environmental  Acoustic  Research  Group  at  the  Naval 
Postgraduate  School  is  engaged  in  research  to  establish  beneficial 
and  detrimental  environmental  effects  important  to  present  and 
future  Navy  acoustic  systems. 

Pursuant  to  the  above  objectives,  environmental  and 
acoustic  models  are  used  to  interpret  and  predict  the  complex 
results  obtained  when  actual  experimental  or  operational  scenarios 
are  utilized.  Acoustic  parabolic  wave-equation  models  are  useful 
for  long-range  environmental  acoustical  studies.  However  the 
basic  assumptions  and  certain  errors  in  the  use  of  these  models 
are  not  always  obvious.  For  this  reason  the  present  work  was 
initiated  to  provide  a  tutorial  introduction  mainly  for  use  of 
students  in  the  Environmental  Acoustic  Research  Group  at  the  Naval 
Postgraduate  School  and  others  interested  in  obtaining  an  initial 
orientation  in  this  field  of  research. 

Presently  five  professors  from  the  Departments  of 
Oceanography  and  Physics  are  involved  in  this  research  as  well  as 
approximately  ten  graduate  students. 


I. 


The  Parabolic  Wave  Equation 
A.  Introduction 


Historically,  the  ( Leontovich-Fock)  parabolic- equation 
approximation  was  developed  for  dealing  with  electromagnetic 
propagation;  see,  for  example,  Fock  (1965).  This  approximation  to 
the  elliptical  wave  equation  made  its  way  into  the  area  of 
acoustic  propagation  through  the  work  of  Hardin  and  Tappert 
(1973),  who  applied  the  split-step  Fourier  transform  method  for 
its  evaluation. 

Since  its  introduction,  the  parabolic  equation  has  'ieea 
subjected  to  substantial  analysis  in  the  acoustics  community,  and 
a  variety  of  alternative  methods  for  solving  the  equation  have 
been  developed.  It  is  the  purpose  of  this  report  to  present  a 
fundamental  introduction  to  the  ’-arabolic- equation  apr roximati^n 
with  some  discussion  of  the  more  viable  methods  of  numerical 
solution.  Some  of  the  advantages  and  disadvantages  of  these 
methods  will  be  noted,  but  detailed  discussion  of  the  problems  and 
requirements  in  implementing  these  methods  ac  computer  algorithms 
will  not  be  treated;  these  lie  beyond  the  simple  introduction 
attempted  herein. 

For  an  overview  of  these  and  other  methods  of  solving  the 
wave  equation  and  its  various  approximate  forms,  a  good  starting 
point  is  DeSanto  (1979),  and  in  particular  the  articles  contained 
therein  by  DeSanto,  and  DiNapoli  and  Deavenport. 
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B. 


Derivation 


If  an  acoustic  pressure  field  of  constant,  angular  frequency 
u)  =  27Tf  is  assumed,  then  the  source-free  linear  acoustic  wave 
equation 


V2E  =  C_2(^£/^t2) 


where  the  speed  of  sound  c  is  a  function  of  space  reduces  to  the 
Helmholtz  equation  (often  termed  the  "time  independent"  wave 
equation) .  For  the  spatial  factor  S  of  the  pressure 

£  =  S(space)  exp(-iaJt)  1.2a 


this  becomes 


where 


72s  +  k2S  =  0 


1.2b 


k  =  4U/c(space)  1.2c 

If  we  now  assume  cylindrical  symmetry  and  adopt  cylindrical 
coordinates  (r,z)  where  z  represents  depth  and  r  is  the  horizontal 
distance  from  the  z-axis,  (1.2)  becomes 


d2S/dr2  +  r^S/di)  +  £ 2S/d  z 2  +  k2S  =  0 


It  is  useful  to  define  k  in  terms  of  an  arbitrary  constant  value 
k0  and  the  index  of  refraction 

n  =  c0/c(r,z)  where  w/k0=Co  1.4 

so  that. 

k2  =  n2kR 


1.5 


Let.  us  now  write  the  solution  for  S  in  the  form 


S  =  H^1]  (k0r)  f  (r,z) 


1.6 


which  when  substituted  into  (1.3)  results  in  the  equation 
t2f/dr2  +  d2f/t)z2  + 


1  +  2 
r  „(l) 


dH^1} (k  r) 
o  o 


H'  (k  r) 
o  o 


dr 


(if  +  (k2-k2)f  =  0 
— «■  o  — 

■jV 


1.7 


If  we  now  restrict  attention  to  ranges  r  such  that 


Qr  >>  1  and  (kQr)  ft#  \fl/  (Trk0r )  exp[  i  (kQr- 77/4  ] 


1.8 


then  [  ]  simplifies  to  2ikQ  and  (1.7)  becomes  well-approximated  by 


J2i  2  2 

+  - — +  2ik  —  +  (k  -k  )  f  =  o 
2  °^r 


c)  z2  ^ 


1.9 


The  pivotal  assumption  that,  f  varies  slowly  with  respect,  to  range, 


J  J  <<:  |2k0^f/<^r 

results  in  the  "parabolic  (wave)  equation" 
a2f/dz2  +  2ik0Of/^r)  +  (k2-k2)f  =  0 


The  assumption  (1.10)  has  some  extremely  important  implications: 

(a)  The  Helmholtz  equation,  an  elliptic  equation,  has  been 
reduced  to  a  parabolic  equation.  This  means  that  the  entire 
acoustic  field  need  not.  be  solved  for  all  relevant  ranges  and 


depths  "simultaneously"  subject  to  boundary  conditions  on  a 
surface  surrounding  the  volume  of  interest.  Instead,  for  an 
outwardly  progressing  wave,  an  initial  boundary  condition  can  be 
established  for  some  small  r  and  then  solutions  for  larger  r's 
obtained  by  increasing  r  incrementally.  This  offers  the 
possibility  of  a  substantial  saving  of  computer  time  and  memory. 
However ,  the  boundary  condition  (initial  condition)  assumed  for 
the  first  range  must  be  carefully  selected.  This  will  be 
discussed  later. 

(b)  The  parabolic  approximation  is  equivalent  to 
neglecting  any  back-scat.t.ering,  since  the  solution  at.  some  range 
r-^  is  the  source  for  the  solution  at  some  larger  range  r2  and 

is  independent  of  any  intervening  changes  in  the  speed  of  sound, 
depth  of  the  water  column,  etc:  As  the  solution  is  stepped  out, 
changes  at  larger  r  can  have  no  effect  on  the  fields  previously 
obtained  tor  smaller  r. 

(c)  The  parabolic  assumption  can  introduce  unavoidable 
errors  in  the  details  of  the  resultant  acoustic  field.  We  shall 
demonstrate  this  by  considering  a  particularly  simple  acoustic 
model  which  can  be  solved  by  both  the  Helmholtz  and  parabolic 
equations.  Comparison  of  the  respective  solutions  will  aid  in 
revealing  some  of  the  inherent  errors  resulting  from  the 
parabolic-equation  approximation . 

For  use  later,  notice  that  for  large  kQr  (1.6)  becomes 

S  =  V2/(Tk0r)  f(r,z)  exp[  i  (k0r-TT/4 )  ]  1.12a 
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and  that  the  acoustic  pressure  in  the  parabolic-equation 
approximation  has  the  form 


£  =  A  V2/ (Trk0r)  £  (  r  ,  z )  exp[ i (kQr-~/4 ) ]  exp(-i'^t) 


1.12b 


with  A  an  arbitrary  constant  and  £{r,z)  the  solution  of  (1.11). 
Recall  this  is  valid  for  kQr  >>  1. 

C.  The  Range-Independent.  Case 

Let  us  assume  that  the  speed  of  sound  is  a  function  only  of 
depth  and  that  any  relevant  boundaries  arc  also  range  independent. 
Direct  solution  of  the  Helmholtz  equation  (1.2)  results  in  the 
well-known  summation  of  normal  modes 

(1) 


£m=  exP  ti'Jt)  I  A*  Vz)  Ho  (kmr) 
m 


1.13a 


where  the  constants  ^  are  determined  by  the  properties  of  the 
Zm  and  source  depth,  and  the  depth-dependent,  functions  Zm  are 
solutions  of  the  equation 


d2Z 


m  + 


dz 


[it2  (*>-**] 


z  =  0 
m 


1.13b 


The  eigenfunctions  Zm  and  the  eigenvalues  km  are  established 
by  the  function  k(z),  the  boundary  conditions  at  the  top  and 
bottom  of  the  water  column,  and  the  properties  of  any  ocean  bottom 
(if  important).  In  the  limit  of  large  r  each  normal  mode  has 
asymptotic  behavior 


^-►expHwt^Zjz)  V/2/(TTkmrj,exp[’  Kk^r-  nj 

1.  ’  4 


1.14 
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Solution  of  the  same  problem  for  the  parabolic  wave  equation 


proceeds  analogously:  Let. 


f  =  R  (r)  Z  (z) 
— m  — m  m 


1.15 


Substitute  into  (1.11)  and  perform  the  usual  separation  of 
variables.  The  result  is  the  pair  of  equations 


d2Z  /dz2  +  [k2( z) -k2]Z  =  0 

m  mm 


1.16a 


and 


2ik  (dR  /dr)  +  (k2-k2)R  =  0 

o  —  m  m  o  — m 


1.16b 


Note  that  (1.16a)  and  (1.13b)  are  identical. 

Equation  (1.16b)  can  be  solved  by  direct  integration  to 

yield 


=  exP 


/  *2-*2  \ 
(  1  —2 — —  r  ) 

\  2k  J 

'  Of 


1.17 


Collection  of  terms  reveals  the  parabolic-equation  solution 

m  /  k2-k^  \ 

P*  =  exp(-iut)  T  \  Zm<z)  Ho  (kor)  exP  (  i'~2L-2~  r 
m  \  2k 


In  the  limit  of  large  r,  each  term  has  asymptotic  form 

_Zz  r  ,2.,  2  “j 


1.18 


P ' 
— m 


exptotjA,,  i(-| T^r-  V/A) 

l.  o  J  _  o 


1 1 

J 


1.19 


( 


N 


:| 

l 


From  (1.14),  we  see  that,  the  phase  speed  cm  of  the  m-th  normal 
mode  is 


c  =  CJ/k 
m  m 


1.20 


whereas  from  (1.19),  the  equivalent,  term  solving  the  parabolic 
approximation  is 

°;  -  u2ko  /  +  ko»  i-2i 

Were  the  acoustic  field  to  consist,  of  only  a  single  mode, 
m  =  M,  then  the  choice  kQ  =  kM  would  exactly  eliminate  the 
phase  speed  error,  and  the  solution  of  the  parabolic  equation 
would  be  identical  with  that,  of  the  Helmholtz  equation. 
Unfortunately,  this  is  not.  usually  the  case.  However,  since  kQ 
is  an  arbitrary  constant,  it.  is  clear  that  if  the  acoustic  field 
is  composed  of  a  set.  of  normal  modes  whose  values  of  km  all  lie 
very  close  together  (a  "narrow  band"  of  modes),  all  other  modes 
being  negligibly  small  or  absent,  then  the  choice  kQ  =  <km> 
where  the  average  <>  is  taken  over  just,  this  narrow  band  will  tend 
to  minimize  phase  errors.  Even  here,  however,  the  errors  may  not. 
be  trivial. 

First,  we  see  that  except  for  the  special  case  k0  =  kM 
the  phase  speed  for  the  normal  mode  is  different,  from  the 

analogous  phase  speed  c^,  for  the  equivalent,  term  in  the 
parabolic  solution.  This  means  that  the  spatial  pattern  of  the 
phase -coherent,  combination  of  the  pressure  terms  will  be 
distorted.  There  is  an  additional  and  equally  important  effect, 
resulting  from  the  fact,  that  the  phase  speeds  for  individual  modes 
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do  not  change  proportionately  with  those  for  the  equivalent  terms 
in  the  parabolic  equation.  Differentiation  of  cm  with  respect 
to  km  yields 


1.22 


1.23 


Thus,  the  analogous  phase  speeds  are  not  related  simply  by  a 
constant,  which  would  merely  "stretch"  or  "shrink"  the  entire 
interference  pattern  with  respect  to  range.  Instead,  the  detailed 
interference  pattern  of  the  acoustic  field  will  also  be  changed. 
For  a  clear  example  of  these  effects,  see  Fig.  5  of  McDaniel 
(1975-1) . 

An  estimate  of  the  maximum  range  for  which  the  parabolic 
equation  retains  sufficient  accuracy  can  be  obtained  following  the 
development  of  Fitzgerald  (1975): 

Assume  that  the  acoustic  field  is  made  up  primarily  of  a  set 
of  strong  modes  with  indices  m  lying  between  m(max)  and  m(min). 
This  could  correspond  to  a  field  consisting  of  trapped  modes  in  a 
mixed  layer,  a  shallow-water  channel  with  a  fast  bottom,  or  the 
deep  sound  channel. 

To  minimize  error,  set 

^o  =  <^m> 


or,  almost  equivalently 


since  either  will  be  about  half  way  between  the  values  found  for 
n(max)  and  m(min) .  Define 


and 


Mc{ 


'm(raax) 


A  kV 


m(max) 


«  k 


m  (min) 


1.24a 


1.24b 


The  accuracy  requirement  can  be  approximated  by  the  condition 


m _ o 

2k 


r  -  k  r 
m 


2 


in  -  m  (max)  or  m(min)  1.25a 


Elementary  manipulation  with  (1.24)  yields 

r  <  TTk0/(Ak)2 


1.25b 


or,  with  the  help  of  c£c0«*c3, 

r  <  Trc^/tuXAc*  )  2]  .  1.25c 

Equation  (1.25)  is  somewhat  less  restrictive  than  Fitzgerald’s 

result,  but  serves  as  a  reasonable  guideline.  Note  the  explicit 
frequency  dependence  in  (1.25c):  For  a  given  family  of  excited 
modes,  the  maximum  permissible  range  for  given  accuracy  will 
decrease  with  increasing  frequency. 

D.  Improving  the  Accuracy  of  the  Parabolic  Equation 
(1)  The  "Pseudoproblem" 

In  the  light  of  the  phase-coherence  difficulties 
discussed  above.  Brock  et  al.  (1977)  investigated  the  feasibility 
of  modifying  the  problem  to  reduce  these  effects.  They  were 
guided  by  ray-tracing  predictions  of  the  turning  points  of  the 


rays  and  the  requirement  of  matching  the  normal  mode  and  parabolic 
phase  speeds,  at  least  over  the  narrow  band  of  1'm's  required  for 
the  validity  of  the  parabolic  equation  (l.ll).  Based  on  these 
considerations,  they  determined  that  an  approximate  analogous 
"pseudoproblem"  could  be  formulated  for  which  the  index  of 
refraction  and  the  depth  at  which  a  particular  speed  of  sound  was 
found  were  adjusted  according  to  the  mapping 


(n,z)-*-(n*,z*) 

where 

ij 

n*  ft?  ( 2n  -  1) 

ij 

z*jfc-  n  z 


1.26a 

1.26b 

1.26c 


Utilizing  this  technique,  they  determined  that  sensitivity  to  the 
choice  of  cQ  was  considerably  reduced,  and  the  predictions  of 
the  transmission  loss  given  by  the  "pseudoproblem"  matched  those 
predicted  by  the  normal  mode  solution  of  the  original  problem  much 
better  than  the  solution  of  the  parabolic  equation  without  the 
mapping  of  (n,z). 

While  these  results  were  obtained  for  the  range- 
independent  case,  the  authors  make  the  plausible  assertion  that  if 
in  a  range-dependent  problem  the  mapping  is  to  be  done  at  each  new 
range  step,  then  the  improvement  of  results  over  those  without  the 
mapping  should  be  about  the  same  as  for  the  range- independent 
case. 

(2)  Alternative  Equations 

It  should  be  pointed  out  that  (1.11)  is  not  the  only  form 


that  a  parabolic  approximation  to  the  wave  equation  can  take. 


Indeed,  several  investigators  have  made  rather  extensive  studies 
of  alternative  forms  and  more  accurate  approximations.  We  shall 
confine  our  discussion  but  provide  references  for  the  reader 
interested  in  pursuing  these  extensions  further. 

McDaniel  (1975—2)  studied  several  methods  of  separating 
solutions  to  an  asymptotic  wave  equation  into  outgoing  and 
incoming  components.  Depending  on  t’  e  method,  when  any 
back-reflected  (incoming)  component  is  neglected,  modified 
parbolic  equations  result  which,  when  compared  to  the  asymptotic 
wave  equation,  reveal  errors  of  various  orders.  Of  the  three 
cases  studied,  two  led  to  second-order  errors  and  one  led  to 
fourth-  order  errors .  The  commonly-encountered  parabolic  equation 
(1.11)  was  one  of  the  second-order  cases.  In  addition,  numerical 
analyses  using  different  algorithms  were  performed  and  results 
were  checked  for  internal  consistency  by  exchanging  source  and 
receiver  positions  and  verifying  that  acoustic  reciprocity  held  to 
reasonable  accuracy. 

Palmer  (1976)  investigated  improvements  to  the  '  ikonal 
equation  and  approximations  to  normal-mode  theory  by  assuming  that 
the  Fikonal  equation  can  be  applied  in  the  horizontal  pla^e.  This 
leads  to  expressions  for  the  normal  mode  coefficients  and  the 
development  of  an  appropriate  Green's  function.  After  rather 
elaborate  mathematical  development,  some  modified  parabolic 
equations  can  be  extracted.  The  t1  rust  of  the  discussion, 
however,  is  toward  a  further  understanding  of  the  plausibility  of 
the  physical  restrictions  necessary  to  justify  the  validity  of 
discarding  small-order  terms  in  the  Helmholtz  equation  to  obtain 


a  parabolic  equation. 


An  investigation  by  DeSant.o  (19  77)  into  the  mathematical 
relationship  between  the  solution  to  the  Helmholtz  equation,  S, 
and  the  solution  _f  to  the  parabolic  equation  (1.11)  yielded  a 
collection  of  correction  terms.  DeSant.o' s  approach  was  to  assume 
that  £5  and  could  be  related  by  an  integral, 

r*° 

S(r,z)  =  AQ  \  f (y,z)  R(y,r , z)  exp [B (y , r) ] dy  1.27 


o 

where  Aq  is  constant,  R  and  B  are  unknown  functions,  and  y  is 
the  dummy  variable  of  integration.  (The  above  integral  is  for 
cylindrical  coordinates,  a  special  case  of  the  more  general 
formulation  accomplished  by  DeSant.o.)  If  this  is  substituted  into 
the  Helmholtz  equation  and  k(r,z)  written  in  the  form 


k(r,z)  =  k1(z)  +  k2(r,z)  where  »  k2 


1.28 


then  it  is  possible  to  obtain  B  and  the  functional  dependence  of  f 
on  y  by  requiring  self-consistency.  What  remains  is  a 
differential  equation  for  R.  DeSant.o  then  shows  that  the  solution 
f?  to  the  parabolic  equation  (1.11)  results  from  the  stationary- 
phase  approximation  of  the  integral  (1.27).  Retaining  higher 
accuracy  in  evaluating  the  integral  provides  correction  terms  to 
the  parabolic  equation  and  therefore  to  f.  In  a  later  paper, 
DeSant.o,  Perkins,  and  Baer  (1978)  begin  with  a  "corrected 
parabolic  approximation"  [compare  with  (1.12a)] 

S  =\/2/ (TTkorV  exp  (kQr-7r/4)] j*f+ ( ir/ 2k)  (<)2f /<)  r 2)J  1.29 


II.  Boundary  Conditions  for  the  Parabolic  Equation 

Except  for  special  cases,  the  parabolic  equation,  like  the 
Helmholtz  equation,  ^oes  not  yield  analytical  solutions  for  space- 
dependent  speed  of  sound  profiles  or  irregular  boundaries. 

Instead,  numerical  methods  must  be  adopted.  There  are  r.  number  of 
numerical  techniques  now  available  for  use  with  the  parabolic  wave 
equation.  We  shall  mention  some  of  f'ose  fat  are  currently 
popular. 

The  major  drivinq  forces  developing  computer  algorithms  for 
numerical  solutions  are  that  computers  are  limited  in  available 
memory  and  computer  time  is  expensive.  As  a  result,  emphasis  has 
been  placed  on  fast-running  programs  which  require  relatively 
little  memory.  Since  the  rarabolic  equ--t:on  is  designed  to  be 
stepped  out  in  range,  it  is  important  to  use  techniques  which 
allow  the  largest  possible  increments  in  both  depth  and  range. 
Since  each  step  requires  numerical  mathematical  manipulation  of 
input  data  and  the  results  of  the  previous  range  step,  efficient, 
computational  schemes  are  required.  In  this  report  we  say  little 
about  these  aspects  of  the  problem;  our  purpose  is  to  describe  the 
methods  rather  than  discuss  the  details  of  their  advantages  or 
disadvantages  as  far  as  computer  implementation  is  concerned. 

Before  turning  to  the  models,  it  is  necessary  to  discuss  two 
aspects  of  the  parabolic  equation  which  are  common  to  all  methods 
of  solution.  In  every  case,  it  is  necessary  to  begin  the 
computation  with  an  input  data  set  of  the  values  of  £_  at  some 
initial  range  as  a  function  of  depth.  This  is  the  initialization 
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problem.  The  second  aspect  is  that,  of  treating  the  boundaries  of 
the  water  column  and  bottom.  There  is  little  difficulty  with  the 
surface,  which  is  represented  by  a  pressure  release  boundary  so 
that,  at  zero  depth  _f(r,0)  =  0  for  all  r.  However,  the  bottom,  if 
an  important,  aspect  of  resultant,  acoustic  field,  presents  some 
difficulties . 


A.  The  Initialization  Problem 

Whatever  set.  of  values  for  _f(rQ,z)  at  the  initial  range 
rQ  are  chosen,  they  must,  be  consistent,  with  the  acoustic  source 
generating  the  pressure  field.  It.  is  therefore  useful  to  obtain  a 
few  results  for  an  omnidirectional  point,  source  located  at  (0,Z) 
in  an  ocean  whose  properties  vary  only  with  z.  This  "boundary 
condition"  can  be  built,  into  the  acoustic  wave  equation  by 
including  a  "source  term". 

If  the  omnidirectional  source  has  unit,  pressure  amplitude  at. 
a  distance  of  1  m,  then  the  appropriate  inhomogeneous  wave 
equation  in  cylindrical  coordinates  is 

[v2  +  k2  (z)  1  S  =  -  2$(r)  $(z-Z)/r  2.1 


The  presence  of  the  term  on  the  RHS  guarantees  that. 

lim  S  =  [r2  +  (z-Z)2]  **  exp[ikVr~2  +  (  z-Z) 2  ] 

( r , z)-K  0  ,  Z) 


2.2 


Given  that,  conditions  exist  for  the  trapping  of  sound 


in  a  channel,  it.  is  plausible  t.o  perform  the  separation  of 
variables 


Vr)zm(z>  2-3 

m 

and  assume  that,  t.he  Zm  satisfy 

d2Zm/dz2  +  [cu2/c2  (z)  -  k2]Zm  =  0  2.4 

the  appropriate  boundary  conditions , and  are  normalized.  Then, 
the  Zm  form  an  ort.honormal  set.  of  eigenfunctions.  Substitution 

of  these  results  into  (2.1)  yields 

H  ["  d2R  /  dr2+r_1( dR  /  dr)+k2R  1  Z  =  -2$(  r)  S  (  z-Z) /r 

m  I-  — m  — m  m— mj  m 

2.5 

If  both  sides  of  this  equation  are  multiplied  by  Zn(z), 
integrated  over  z,  and  use  made  of  t.he  ort.honormalit.y  condition 

[ Z  Z  dz  =  b  26 

I  m  n  mn  • 0 

then  t.he  result,  is 

d2R  /dr2  +  r-1(  dR  /dr)  +  k2R  =  -  2  <5(  r )  Z  (Z)/r  2-7 

— m  — m  m— m  m 

which  is  solved  by 

R  =  ilT Z  ( Z )  H  ( 1)  (k_r)  2.3 

~ m  mom 

Now,  if  (2.7)  is  substituted  back  into  (2.5),  we  obtain  the 
useful  relationship 


This  expression  ignores  the  collection  of  continuous 
eigenfunctions,  since  these,  practically  speaking,  contribute 
nothing  to  the  acoustic  field  at  ranges  of  interest. 

(1)  The  Gaussian  Field. 

When  an  omnidirectional  source  of  sound  is  reasonably 
distant  from  either  the  ocean  surface  or  any  bottom,  then  the 
source  may  be  approximated  by  a  Gaussian  pressure  field  for  the 
initial  set  of  pressure  as  a  function  of  depth.  The  difficulty  is 
to  determine  the  parameters  of  the  Gaussian  distribution  to 
"match"  the  point  source  to  the  particular  propagation  problem 
under  consideration. 

In  the  immediate  vicinity  of  an  omnidirectional  source 
which  is  not  too  close  to  any  reflective  boundary,  the  amplitude 
of  the  pressure  must  decrease  with  distance  from  the  source, 
according  to  spherical  spreading.  For  such  a  source  located  at 
(r,z)  =  (0,Z)  therefore,  from  (1.2)  the  amplitude  of  the  radiated 
pressure  must  be  given  by 

P  =  [r2  +  (z-Z  )2yh  2.10 

(Recall  we  are  in  cylindrical  coordinates  with  radial  symmetry  and 
have  assumed  unit  pressure  amplitude  at  a  distance  of  1  m.  This 
amplitude  choice  facilitates  conversion  between  source  level  SL 
and  pressure  amplitude  for  sources  of  arbitrary  strength.) 

In  the  case  of  an  infinite,  homogeneous  medium  the  source 
with  amplitude  given  by  (2.10)  and  angular  frequency  cj  must  be 
described  by  an  appropriate  collection  of  delta  functions  at  (0,Z) 
However,  all  the  energy  radiated  from  the  source  does  not  find  its 
way  into  the  sound  field  trapped  by  the  sound-speed  profile. 
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Instead,  energy  can  be  lost,  t.hrough  interaction  with  the  bottom. 
Thus,  rays  whose  angles  of  elevation  or  depression  exceed  certain 
limits  will  be  lost  to  the  channel  at.  large  ranges,  and  the  source 
appears  as  if  it.  is  not.  a  point,  source,  but.  instead  possesses 
vertical  directivity. 

Brock  (1978)  presents  one  way  of  incorporating  this 
apparent,  directivity  into  the  parabolic  equation  initialization: 

Assume  that  in  the  volume  of  water  surrounding  the  source 
at.  (0,Z)  acoustic  conditions  are  relatively  uniform.  Then  we  can 
treat  k  as  constant  for  z~Z  and  r-0.  If  we  then  take 

k  =  kQ  =  u)/ c  ( Z )  2.11 

where  c(Z)  is  the  speed  of  sound  near  the  source,  the  parabolic 
equation  (1.11)  simplifies  t.o 

b2t /bz2  +  2iko(df/d^)  =  0  2.12 

Within  the  volume  for  which  (2.11)  is  a  reasonable  approximation, 
but  still  under  the  condition  kQr  >>  1,  we  require  that  the 
amplitude  of  the  solution  t.o  (2.12)  be  consistent  with  the 
amplitude  P'  given  by  (2.10).  We  can  then  extrapolate  f_  back  to 
r  =  0  and  obtain  an  extrapolated  equivalent,  boundary  condition 
along  the  z-axis.  (Physically,  this  seemingly-artificial  approach 
is  equivalent,  to  taking  the  far-field  behavior  of  a  directional 
source  and  extrapolating  it.  back  t.o  the  position  of  the  source, 
thereby  ignoring  near-field  effects.) 


This  extrapolation  is  accomplished  by  reversing  the  steps 


of  the  argument..  Postulate  that  the  equivalent,  boundary  condition 
for  f  at  (0 , z)  is  given  by  the  Gaussian  distribution 


f  (0 ,  z) 


F  exp[  —  (  z-Z)  2/C"2  ] 


2.13 


where  F  and  a~  are  to  be  determined.  With  this  as  an  initial 
condition  (and  requiring  that,  f^  vanish  at.  large  r),  according  to 
Brock  the  solution  of  (2.12)  can  be  verified  to  be 


f(r,z)  =  [F/^(r)  ]exp|-(  z-Z)  2/  'fr(r)  ]| 


2 . 14a 


g_(r)  =  1  +  L  2r/(kQ^2) 

Now,  substitute  (2.14)  into  (1.12)  and,  since  we  have  two 
arbitrary  constants,  set.  A  =  (T  kQ/2)  1/2  for  convenience. 
Next,  take  P  =  |  gj, 

P  =  |f(r,z)|  /VT* 

Brock's  result,  is 

P2  =  (F2/r)  (h/ l/l+h2  )  exp'f-2h2(z-Z)2/[  £T2(l+h2)]^ 


2.14b 


2.15 


2.16a 


with 

h  =  j  kQ(r  2/r  2. 16b 

The  next  step  is  to  expand  (P1  )2  and  P2  in  power  series.  For 
(P*  )2,  assume  that,  [(z  -  Z)/r]2  is  small.  For  P2,  assume 
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that  the  exponent,  in  (2.9a)  is  small.  This  leads  to  two  power 
series  in  [(z  -  Z)/r]^#  and  equating  coefficients  of  the  leading 
two  terms  yields  the  results 


c r  =  y?/k 


2 . 1  /a 


f  =  v'vkVo-  =  yr 


2.17U 


Thus,  the  initial  values  of  f  are  given  by 


f  (0  ,  z)  =  ./jP  exp[-k^(  z-Z)  2 / 2 D 


2.18 


for  an  omnidirectional  source  with  unit,  pressure  amplitude  at  1  m . 

Direct,  substitution  of  (2.17a)  into  the  exponent,  of 
(2.16a)  reveals  that,  the  exponent,  is  small  if  [(z  -  Z)/r]2  is 
small.  Thus,  the  physical  implication  of  the  approximations 
leading  t.o  (2.18)  from  (2.16)  is  that,  the  omnidirectional  source 
i3  approximated  by  an  equivalent,  directional  source  whose  acoustic 
axis  lies  in  the  horizontal  plane  and  whose  beamwidtl.  satisfies 
the  requirement,  of  relatively  small  angles  of  elevation  and 
depression.  It.  must  be  noticed  that,  the  Gaussian-f ield  approach 
may  run  into  difficulties  if  the  source  is  too  close  t.o  the  water 
surface  or  the  water-bot.t.om  interface.  In  either  of  these 
situations,  the  Gaussian  field  may  intercept,  the  boundary  before 
it  has  become  negligibly  small.  As  a  rough  criterion,  the  source 
should  be  removed  by  several  cr  from  any  boundary. 


(2)  Normal  Modes 


While  the  Gaussian  field  initialization  has  the  advantage 
of  beginning  with  a  field  which  varies  rather  smoothly  with  depth 
and,  if  fairly  narrow  in  width,  appears  to  be  a  reasonable 
approximation  to  a  point  source  (at  least  to  the  eye),  it  is  not 
universally  accepted  as  being  without  inherent  error.  See,  for 
example.  Wood  and  Papadakis  (1980).  An  approach  based  on  normal 
modes  can  be  conceptualized  as  follows: 

For  ranges  near  the  source,  outward  propagation  can  be 
described  by  assuming  that  the  speed  of  sound  profile  and  bottom 
properties  are  uniform  throughout  the  medium  (this  includes  the 
depth).  The  Helmholtz  equation  (2.4)  can  be  solved  numerically  to 
obtain  the  depth-dependent,  eigenfunctions  Zm(z).  Since  the  Zm 
form  an  orthonormal  set,  they  can  be  used  to  represent  a  point 
source.  Then,  for  a  point  source  at  (0,Z)  of  unit  pressure 
amplitude  at  r  =  1  m,  we  have  from  (2.9)  and  (2.8) 


■,z)  =  i  If  ^ 


S  (r ,  z)  =  i  TT  /  HU)  (k  r)  Z(Z)  Z  (z) 
•  o  m  m  m 
m 


2.19 


(1.6)  allows  us  to  form  an  initial  expression  for  f_(0,z), 

f  (0,z)  =  i  If)  Z  (Z)  Z  ( z) 

—  1 —  m  m 


2.20 


m 


m 


Because  of  the  small  angular-aperture  assumption  implicit  in  the 
parabolic  equation,  and  also  because  in  most  situations  of 
practical  interest  only  modes  corresponding  t.o  rays  with  small 
angles  of  elevation  and  depression  with  respect  to  the  horizontal 
are  trapped,  the  summation  over  m  can  be  restricted  just  to  that 
subset  of  Zm  satisfying  these  conditions. 
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Notice  that,  since 


only  this  subset  is  to  be  retained,  t.he  Helmholtz  equation  must,  be 
solved  only  t.o  obtain  this  subset,  which  results  in  a  substantial 
saving  in  computer  running  time . 


Given  the  initial  expression  (2.20),  t.he  discrete  depth 
values  zm  are  specified,  t.he  associated  f(0,zm)  found,  and  the 
solutions  for  increasing  ranges  stepped  out.  If  practical  aspects 
of  the  problem  preclude  starting  at.  r  =  0,  the  normal  inodes  will 
have  to  be  allowed  t.o  propagate  out  to  the  minimum  range  rQ  at 
which  the  parabolic  equation  can  be  implemented.  In  such  a  case, 
errors  introduced  as  a  result,  of  approximations  necessitated  in 
utilizing  the  normal  mode  problem  must,  be  studied  carefully.  For 
example,  the  bottom  may  have  to  be  unrealistically  smoothed  out.  to 

ro‘ 

For  an  example  utilizing  a  normal-mode  starter  (and 
displaying  some  of  the  difficulties  encountered),  see  Guthrie  and 
Gordon  ( 1977 )  . 

(3)  Ray  Tracing 

Another  approach  to  obtaining  an  initial  set.  of  values 
for  f^  at  some  finite  range  rQ  from  the  source  involves  phase- 
coherent  ray  tracing.  From  a  practical  point,  of  view,  this  has 
appeal,  since  interactions  of  any  ray  with  the  bottom  can  be  dealt, 
with  by  introducing  t.he  appropriate  plane-wave  reflection 
coefficient,  as  long  as  the  source  is  several  wavelengths  away  from 
the  bottom.  The  required  family  of  rays  is  traced  out  to  t.he 
desired  initial  range,  with  all  phase  information  retained 
(including  that,  arising  from  surface  and  bottom  reflections),  and 
then  combined  to  give  t.he  resultant  pressure  distribution  at 


(rQ,z).  From  this,  f(rQ,z)  can  be  obtained  and  the  f ' s  for 
larger  ranges  found  by  stepping  out  with  whatever  algorithm  is 
used  to  solve  the  parabolic  equation.  For  an  example,  see  Guthrie 
and  Gordon  (1977). 


B.  Treatment  of  the  Bottom 

Adequate  representation  of  the  bottom  hns  been  a  difficulty 
with  the  parabolic  wave  equation.  The  seriousness  of  the  problem 
and  how  it  is  dealt  with  depends  on  :-he  algorithm  used  in  stepping 
the  solution  out  in  range. 

We  can  say  that  many  of  the  approaches  fall  into  three 
categories: 

(a)  Reflection  of  the  equivalent  rays  from  the  bottom  allows 
the  reflective  loss  to  be  given  as  a  function  of  the  apparent 
angle  of  incidence  of  the  ray  on  the  l  ot  .oiu.  This  can  be  either 
specified  as  input,  or  calculated  according  to  some  formula.  The 
possibility  of  the  bottom  depth  >  aing  a  function  of  range  is 
included.  (If  attention  is  restricted  to  sufficiently  large 
range,  of  course,  all  rays  less  grazing  than  critical  can  be 
neglected  since  these  correspond  to  bottom  bounce  paths  with 
appreciable  transmission  into  the  bottom  at  each  interception. ) 
However,  the  determination  of  what  particular  rays  are  striking 
the  bottom  at  each  specified  range  is  a  nontrivial  process,  and 
the  models  and  techniques  used  in  this  determination  often  contain 
rather  sweeping  approximations  and  may  impose  auesti onable 
requirements  on  the  depth  dependence  of  the  speed  of  sound  and  the 
absorption  coefficient  in  the  bottom.  An  example  of  this  approve1'' 
can  be  found  in  the  report  by  Stieglitz  et  al.  (1979). 

(b)  The  speed  of  sound  and  the  absorption  coefficient  can  be 
specified  as  functions  of  depth  by  a  complex  index  of  refraction 
for  the  bottom,  and  this  information  built  into  the  algorithm. 
Depending  on  the  computational  scheme,  it  may  be  necessary  to 


prohibit  any  discontinuity  in  c(z)  at  the  water-bottom  interface 
by  allowing  c(z)  to  change  rapidly  but  smoothly  over  some  small 
interval  in  depth  about  the  interface,  and  to  have  the  bottom 
absorption  coefficient  rise  with  depth  from  zero  at  the  interface 
to  finite  value  at  larger  depths .  In  this  method  it.  may  be 
necessary  to  introduce  an  artificial  pressure-release  boundary  at 
some  appreciable  depth  in  the  bottom.  This  step  is  not 
particularly  troublesome,  if  done  correctly:  For  an  absorbing  or 
non-absorbing  bottom,  the  acoustic  field  will  eventually  begin  to 
decay  quasi-exponent.ially  so  that  the  energy  found  at  all  but 
shallow  depths  in  the  bottom  becomes  quite  small;  if  this  is 
allowed  to  be  reflected  back  up  through  the  bottom  and  into  the 
water,  the  error  introduced  into  the  resultant  acoustic  field  is 
negligible  if  the  reflecting  surface  is  deep  enough.  See,  for 
example,  Williams  (1975)  and  the  references  he  cites. 

(c)  If  the  algorithm  permits,  the  water-bottom  interface  can 
be  built  into  the  computer  code  directly  in  terms  of  the  boundary 
conditions  of  continuity  of  pressure  and  continuity  of  the  normal 
component  of  the  particle  velocity.  Let  the  water  be  labeled 
fluid  1  and  the  bottom  fluid  2,  and  assume  a  flat  bottom  at  depth 
z  =  ZB.  Then  (for  a  flat  horizontal  bottom)  the  boundary 
conditions  on  f?  become 

f](ZB)  =  f2(ZB)  2.21 

and 

=,V1(^2/^)Zb  2.22a 

D  D 
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Notice  that,  if  the  bottom  is  not  flat,  so  that  ZB  is  itself  a 
function  of  range,  (2.21)  remains  the  same  but  the  continuity  of 
normal  particle  velocity  takes  the  more  general  form 


n*  (  Vfx) 


n*  ( 


B 


2.22b 


where  n  is  the  unit,  normal  to  the  bottom  at  coordinate  (r,ZB). 


III.  Methods  of  Solution 

There  are,  of  course  many  different  ways  in  which  solutions 
to  the  parabolic  equation  (1.11)  can  be  attempted.  As  with  the 
Helmholtz  equation,  closed-form  or  analytical  solutions  are 
available  only  for  certain  special  cases.  As  a  consequence, 
solutions  are  accomplished  numerically  with  algorithms  resulting 
in  efficient  computer  use.  As  might  be  expected,  the  more 
complicated  the  case,  the  more  difficult  the  solution.  Certain 
methods  of  solution  are  useful  and  of  sufficient  accuracy  only 
when  the  bottom  is  not  important  to  the  problem,  as  in  the  long- 
range  propagation  in  the  deep  sound  channel  when  the  bottom  of  the 
channel  lies  well  above  the  ocean  floor,  or  when  the  bottom  is 
sufficiently  well-matched  that  it  can  be  treated  by  fairly  simple 
approximations.  More  recent  methods  have  been  developed  which 
allow  more  realistic  inclusion  of  a  bottom.  Our  approach  here 
will  be  to  outline  a  few  of  the  methods  in  roughly  chronological 
order  of  appearance,  but  leaving  detailed  discussion  and  the 
consideration  of  fine  points  and  subtleties  for  the  interested 
reader  to  discover  from  the  references. 

Before  discussing  the  methods,  we  shall  develope  a  single 
formalism  to  unify  the  discussion.  Rearrange  (1.11)  to  isolate 
the  derivative  with  respect  to  range. 


^f/dr  =  !">ii(k2-k2)/k  +  (>ii/k  )  (<)2/<}z2)!  f 


3.1 


Following  the  fundamental  lemma  of  calculus 


Ja  =  lim  -(r+*r,2)  ~  -(r>2) 

3  r  ’  av-*o  A  r 


3.2 


we  can  turn  this  into  an  equation  involving  the  increment  br  in 
range,  which  provides  the  basis  for  the  numerical  incrementation 
of  range. 


f(r+br,z)  =  ll  +  (^r)  •sik^n2-!)  +  )  ( >si/ko )  ( *  2/£  z2 )  f(r,z) 


If  we  define  the  operators 


a  =  >si(k2-k2)/ko 


b  =  i/(2kQ) 
D2  =  c)  2  /  o  z  2 


*siko(n  -1) 


3.4a 

3.4b 
3 . 4c 


then  (3.3)  assumes  the  form 


f(r+Ar,z)  =  I  l  +  (br)a+(£r)bD2~|f  (r  ,z) 


3.5 


The  operator  a  describes  the  refractive  properties  of  the  medium 
(and  may  include  those  of  the  bottom).  It  is  thus  a  function  of 
range,  depth,  and  the  constant  kQ.  The  product  of  operators 
bD2  describes  the  depth  dependence  of  f_. 

If  we  make  use  of  the  expansion  of  an  exponential. 


exp(x)  =  1  +  x  +  0(x  ) 


3.6 
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then  it  is  clear  that,  we  can  write  (3.5)  symbolically  in  the  form 


f(r+^r,z)  =  e*r(a+bD  )f(r#z) 


3,7 


if  I  £  r[a  +  bD2  ]fj  <  <  j  f_| .  The  reason  for  this  representation  is 
that.  it.  will  allow  manipulation  of  the  operators  to  be  brief  and 
simplified.  For  example,  note  that,  the  manipulation 


A  +  B  A  B 
e  =  e  e 


is  identical  with 


1  +  A  +  B  =  (1  +  A)  (1  +  B) 


through  terms  of  first,  order  in  A  and  B.  The  difference  between 
the  two  sides  of  the  above  equation  is  2AB,  which  is  of  second 
order.  As  can  be  seen,  for  small  Ar  t.he  operators  a  and  bD2 
acting  on  J:  produce  terms  which  are  substantially  smaller  then  Jf 
itself,  so  this  symbolic  formalism  is  useful. 

A.  the  Split-step  Fourier  Transform 

This  approach  makes  use  of  the  Fourier  transform  and  its 
inverse  to  eliminate  t.he  second  partial  with  respect,  to  z.  There 
are  two  forms  in  which  this  can  be  done;  we  first,  consider  t.he 
form  which  appears  to  be  used  most,  widely. 

Let.  us  define  the  symmetric  complex  Fourier  transform  pair 

. 

1  C7 

3 .  ba 


ps<  )  - 


=  vnr  /_ 


,  ,  isz 

(  )  e  dz 


f‘1(  )  =  J, 


•O 

•6 


T2K 


,  .  -isz , 

(  )  e  ds 


3 .  bb 


where  (  )  represents  any  reasonably  well-behaved  function  of  s  in 

the  transform  (3.8a)  and  any  reasonably  well-behaved  function  of  z 
in  the  inverse  transform  (3.8b).  (For  the  required  mathematical 
properties  defining  "reasonably  well-behaved",  consult  any 
standard  text  on  transform  theory. ) 

Let  us  apply  the  transform  pair  to  the  RHS  of  (3.7)  after 
utilizing  the  approximation 

e*r(i  +  bD2)  Ur)ae(4r)bD2  j-9 


Now, 


1 

V  2  TT~ 


go  2 

£r  bD  _  isz, 
e  —  f  e  dz 


3.10 


arbD  ,  isz 
e  —  fe 


(1  +  Arb  —  )  f  eisz 
“dz2  ~ 


(1  -  Arbs2)  f  elsz 


*  e  -Ar^s  f  eisZ 


3.11 


Because  the  integration  in  (3.10)  is  over  z,  exp(-Arbs  )  can  be 


factored  out. 


Fs  |exp(ArbD2)  exp  (-Arbs2)  Fg(  f ) 


3.12 


and  the  inverse  transform  (3.8b)  applied, 
exp(  ArbD2)  f  ^  jexp( -Arbs2 )  Fg(  f  )j 


3.13 


It  then  follows  at.  once  from  (3.9)  that 

eUr(W)^expUra)p-ire(-irbsJ)F  (f)] 


— s  — 


3.14 


and  insertion  of  this  into  (3.7)  yields  the  split-step  Fourier 
transform  result 


f(r+/ir,z)  =  exp(hra) F_1 [exp( -irbs2 ) F  (f)] 


3.15 


where  a  and  b  are  given  by  (3.4a)  and  (3.4b).  Notice  that  they 
yield  purely  multiplicative  factors.  Examination  of  (3.15) 
reveals  that  the  algorithm  can  be  viewed  as  dealing  with 
diffractive  effects  by  the  exponent  multiplying  Fg  and 
refractive  effects  separately  by  the  exponent  multiplying  the 
inverse  transform. 

Detailed  investigation  of  the  errors  introduced  by  use  of  the 
exponential  approximations  is  involved  and  results  in  rather 
difficult  expressions  to  interpret.  However,  Jensen  and  Krol 
(1975)  have  investigated  the  results  for  situations  typically 
encountered  in  the  ocean  and  obtained  relatively  simple  criteria: 
For  cases  of  practical  interest,  the  errors  introduced  in 
predicting  the  behavior  of  t_  over  the  range  increment.  Ar  are 
bounded  by  the  larger  of 


J  2(aVcQ)  (dn/dz)Ar  (df/dz )j 


3.16a 
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or 


[(£*;/ cQ)  (dn/d  z)^rj 


3 . 16b 


These  reveal  that  if  t_  is  not  sufficiently  slow  in  its  depth 
variation,  then  the  error  increases  linearly  with  frequency  and 
linearly  with  the  range  increment.  For  sufficiently  slowly 
varying  f,  however,  the  error  depends  quadratically  on  both  these 
quantities. 

When  (3.15)  is  implemented  on  a  computer,  the  Fourier 
transform  is  replaced  by  its  numerical  counterpart,  the  FFT .  This 
results  in  a  fast-running  algorithm  which  deals  with  a  finite 
number  of  terms.  The  FFT  requires  a  finite  depth  over  which  it  is 
applied,  so  an  artificial  bottom  (pressure  release)  must,  be 
inserted  at  some  distance  below  the  water-bottom  interface.  Given 
this  interval  in  z,  the  depth  increment  can  be  selected  for 
application  of  the  FFT.  It.  must  be  noted,  however,  that  running 
time  increases  fairly  dramatically  as  the  number  of  depth 
increments  N  is  increased.  Guthrie  and  Gordon  (1977)  estimate 
that  the  time  will  increase  approximately  as  log  (NN).  This  can 
place  rather  drastic  restrictions  on  the  variation  of  the.  index  of 
refraction  n(r,z)  with  depth,  since  the  number  of  points  in  depth 
must  be  great  enough  to  reproduce  with  fairly  high  accuracy  the 
details  of  the  speed  of  sound  gradients. 

An  additional  problem  with  the  split.-step  Fourier  transform 
approach  is  that,  the  change  in  density  between  water  and  the 
bottom  cannot,  be  considered.  While  the  reflection  coefficient  at 
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at  the  water-bottom  interface  is  considerably  more  sensitive  to 
changes  in  the  speed  of  sound  than  to  changes  in  the  density,  if 
the  change  in  the  latter  is  too  great,  realistic  predictions 
cannot  be  expected  from  the  algorithm.  Attempts  to  surmount  this 
difficulty  by  using  equivalent  rays  and  specifying  the  reflective 
loss  (and  phase  angle)  at  the  interface  have  had  mixed  success. 
Indeed,  pessimism  has  been  expressed  as  to  the  utility  of  this 
algorithm  in  situations  for  which  the  bottom  is  an  important 
element  of  the  propagation  problem?  see,  for  example,  Jensen  and 
Krol  (1975)  and  Brock  (1978).  Further,  inclusion  of  a  bottom 
increases  running  time  considerably.  Jensen  and  Krol  (1975)  have 
compared  computer  times  for  split-step  parabolic,  normal  mode,  and 
ray  programs  for  deep  water  (for  which  the  bottom  is  not 
important)  and  shallow  water  for  which  the  bottom  is  important. 
There  is  considerable  variation  in  running  times,  but  in 
particular  for  their  case  the  shallow  water  problem  led  to 
extremely  long  running  times  for  the  split-step  approach. 

The  second  form  of  the  split-step  Fourier  transform  approach 
has  been  described  by  McDaniel  (1975).  This  is  based  on  an 
alternative  formulation  of  the  RHS  of  (3.7)  presented  by  Tappert 
and  Hardin  (1974), 

exp[£r(a+bD  )]«  exp(ara/2 ) expCbrbD  )exp(-ira/2) 

Application  and  manipulation  of  the  Fourier  transform,  its 
inverse,  and  the  exponential  terms  proceeds  much  as  before,  and 


the  details  are  uninteresting.  The  result  is 

f(r+Ar  ,z)«eAr-/2F"1|e_ilr^sF  'e^r-/2  f  SI  3.17 

Comparison  of  this  result  with  the  analogous  (3.15)  reveals  that 
the  index  of  refraction  has  been  treated  a  little  more  carefully. 
When  (3.17)  is  analyzed  for  the  dependence  of  errors  on  the  range 
increment,  it  appears  that  the  error  depends  on  (£r)3  which  is 
an  improvement  over  the  errors  resulting  from  the  approximations 
leading  to  (3.15).  What  restrictions  must  be  placed  on  !>z, 
however,  by  the  introduction  of  n(r,z)  into  the  Fourier  transform 
over  z  do  not  seem  to  have  been  isolated. 
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B.  Alternatives  to  the  Split-step  Method 
Within  t.he  last  few  years,  there  have  been  a  number  of 
alternative  approaches  to  developing  efficient  algorithms  for 
solving  the  parabolic  equation  (1.11).  We  will  mention  three 
closely  related  articles  which  provide  a  rather  interesting  line 
of  evolution.  All  of  these  discard  the  spiit.-st.ep  technique  with 
its  associated  use  of  the  FFT  direct  and  inverse  numerical  method. 
Instead,  the  second  derivative  with  respect  to  depth  is  calculated 
numerically  from  the  depth  increments  by  the  use  of  a  central 
finite-difference  approximation.  Since  the  approaches  based  on 
this  technique  require  more  detailed  specification  of  the 
coordinates  of  any  spatial  point  on  the  range  and  depth  mesh,  we 
will  adopt,  a  convenient,  formalism  which,  although  somewhat 
unconvent.ial,  is  easily  followed  and  succinct..  Since  t.he  depth 
below  t.he  ocean  surface  can  be  written  as  zm  -  m  £z  where  JSz  is 
t.he  chosen  depth  increment.,  and  t.he  range  can  be  written  as  rn  38 
n  A.r,  let.  us  define  t.he  shorthand 

f(r  ,z  )  =  f(nar,m/3z)  5  f[n,m]  3.18 

—  n  m  — 

(Note  that,  n  need  not.  begin  with  n  =  0.  If  the  parabolic  equation 
is  initialized  at  some  finite  range,  then  t.he  appropriate  finite 
value  of  t.he  initial  n  can  be  used.  )  The  underbar  on  the  RHS  has 
been  suppressed  for  notat.ional  convenience.  It.  is  to  be 
understood  that.  f[n,m]  is  complex. 
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Now,  it.  is  straightforward  to  see  that  the  second  depth 
derivative  can  be  estimated  numerically  by  the  formula 


2>2f  ! 


f [n ,m+l ] 


r  tr 


-  2f[n,m]  +  f[n,m-l] 
(az)  2 


3.19 


By  use  of  (3.19),  Lee  and  Papadakis  (1930)  cast  (1.11)  into 
the  form 


df[n,m]/dr  =  a[n,m]f [n,m]  +  [b/(2 z) 2 ] ^ f [n,m+l]  -  2f[n,m]  +  f[n,m- 

3.20 

where  the  operators  a  and  b  have  been  defined  previously,  (3.4). 

(It  is  worth  recalling,  for  this  discussion  and  what  follows,  that 
b  is  simply  a  constant  and  that,  a  depends  functionally  on  depth 
and  may  also  depend  on  range  through  the  square  of  the  index  of 
refraction.  )  This  equation  reveals  that  the  first,  range 
derivative  of  f  at  some  point.  [n,m]  depends  on  the  values  of  f  at. 
the  adjacent  mesh  points  characterized  by  [n,m+l],  [n,m],  [n,m-l] . 
Thus,  (3.20)  provides  a  set  of  first-order  ordinary  differential 
equations,  each  of  which  is  coupled  t.o  its  immediate  neighbors  in 
depth.  While  this  may  appear  t.o  require  a  considerable  amount,  of 
storage  when  these  matrices  are  programmed  into  the  computer,  the 
fact  that  the  coupling  involves  only  adjacent,  depth  values 
collapses  the  matrix  representing  the  RHS  of  (3.20)  to  a 
"tridiagonal"  form  for  which  all  elements  are  zero  except,  those 
whose  depth  indices  are  m,  m  +  1.  The  depth  incrementation  is 
halted  at  the  bottom,  and  Lee  and  Papadakis  let  the  element.  f[n,M+l] 
represent  the  value  of  f  at.  the  bottom. 
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Note  that  M  may  itself  be  a  function  of  n,  so  that  a  sloping  or 
irregular  ocean-bottom  interface  is  allowed.  The  range-depth  mesh 
does  not  project,  into  the  bottom.  In  this  sense  this  approach  is 
incomplete,  but  the  authors  avoid  this  difficulty  to  some  extent 
by  invoking  a  boundary  condition 

Vfc-n  +  ^(s)£  =  0  3.21 

where  n  is  the  local  normal  to  the  interface  and  s  specifies  the 
line  of  the  interface  in  (r,z)  space.  It  can  be  quickly  seen  for 
monofrequency  sound  that  this  is  equivalent  to  the  condition 

-icou-ni  +  of  (s)  £  =  0  3.22a 

which  can  be  rearranged  to  yield  the  form 

£/(u*n)  =  i«**/c/(s)  3.22b 

Thus,  i&Vs<  represents  the  ratio  of  the  pressure  p  to  the 
component  of  the  particle  velocity  u  locally  normal  to  the  bottom. 
For  the  case  of  a  flat  horizontal  bottom,  (3.21)  takes  the  form 

^£/dz  +  2^(r)£  =  0  at  the  bottom  3.23 

Equation  (3.23)  can  be  solved  by  £  =  exp(-«4z),  and  the 
numerical  equivalent  for  specifying  the  interface  value  f[n,M+l] 
is 

f[n,M+l]  =  f  [n,M]  exp  |  -«[n]Az"|  3.24a 
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At  the  ocean  surface,  of  course,  the  pressure  must  vanish  and  the 
surface  boundary  condition  is 

f [n , 0]  =  0  3.24b 

The  authors  solve  two  cases:  a  per fectly-rigid ,  flat,  horizontal 
bottom  and  a  perfectly-rigid,  flat,  sloping  bottom.  For  both,  = 
Results  are  compared  respectively  with  a  normal  mode  solution 
(horizontal  bottom)  and  a  method  of  images  solution  (sloping 
bottom),  and  are  found  to  be  in  very  good  agreement. 

In  a  later  paper,  Lee,  Botseas,  and  Papadakis  (1981) 
reinvestigated  the  solution  of  (3.20),  this  time  using  an  implicit 
finite  difference  method  rather  than  the  linear  and  nonlinear 
multistep  methods  used  in  the  Lee  and  Papadakis  (1980)  article. 

The  practical  advantage  of  the  implicit  finite-difference 
technique  is  that  (at  the  present  time)  while  it  is,  according  to 
the  authors,  about  equally  as  fast  as  solving  (3.20)  by  the 
multistep  methods,  it  requires  signif icent.ly  less  memory.  Several 
propagation  examples  are  worked  out,  using  range  independent, 
upslope,  and  downslope  cases.  The  bottoms  were  fast  and 
absorptive  with  densities  up  to  50%  greater  than  found  in  the 
water  column.  Results  were  compared  with  a  normal  mode  program 
(SNAP)  and  a  split-step  parabolic  equation  (PAREQ),  both  from 
SACLANT.  Propagation  loss  curves  for  all  three  programs  are  shown 
for  the  various  cases,  but.  detailed  comparisons  and  discussions 
were  not  attempted.  Again,  the  boundary  conditions  were  handled 
according  to  use  of  (3.24a)  for  the  horizontal 


bottom,  that  equation  s  generalization  for  the  range  dependent 
bottoms,  and  the  surface  was  represented  by  (3.24b).  The  authors 
4-lso  point  out  the  time-saving  advantage  of  not  having  to  run  the 
parabolic  solution  into  the  bottom. 

McDaniel  and  Lee  (1982)  attacked  the  boundary  condition  at 
the  bottom  more  directly.  If  the  bottom  is  found  at.  depth  m^-'z, 
assumed  to  be  flat  and  horizontal,  and  to  have  properties 
independent  of  range,  then  application  of  continuity  of  pressure 
and  the  normal  component  of  particle  velocity  across  the  interface 
yields 

f  1  [n  =  f2[n,mb]  3.25a 


-1 


3.25b 


where  the  subscripts  1  and  2  refer  to  the  water  column  and  ocean 
bottom  respectively.  After  performing  some  Taylor  series 
expansions  for  the  second  derivative  of  f  with  depth  in  both  media 
and  making  use  of  the  boundary  conditions,  the  authors  are  able  to 
write  a  parabolic  equation  which  must  be  satisfied  at  the 
interface. 


3.26 


/  Pi '  /  Pi  \ 

dfiTi,mb]{l+  ~p2j=[a l^'V  +  J~2  a2ln,mb]j  u[n'mb] 

C  p  i  /  pi ' 

+  —  -*■  ■<  — —  u  [n,m, +1]  -  (  1+  —  u  [n,m.  ]  +  u  [n,m. -1) 

(K*)2  IP  2  b  \  P  2  13  b 


dr 


Away  from  the  interface,  of  course,  (3.20)  nust  hold  in  both 
media.  Thus,  (3.26)  is  solved  for  m  =  i  ^  but  (3.2U)  is  solved 
for  all  other  m.  With  the  bottom  included  in  the  problem,  througn 
application  of  the  boundary  conditions  (3.25),  it  i''  necessary  to 
allow  the  depth  z  to  penetrate  at  least  one  increment  A  z  below  the 
interface.  With  the  single  exception  of  (3.26)  for  m  =  n^,  the 
problem  is  run  according  to  any  of  the  above-mentioned  methods;  if 
the  split-step  FFT  method  is  used,  however,  it  must  be  assumed 
that  there  is  no  density  change  from  water  into  bottom.  An 
example  is  worked  out;  a  two  gradient  concave  speed  of  sound 
profile  overlies  an  isospeed  fast  bottom  with  density  slightly 
more  .than  twice  that  of  the  water.  The  problem  is  solved  using  a 
normal  mode  program,  the  split-step  FFT  method,  and  implicit 
finite-difference  methods  with  and  without  the  interface  condition 
(3.26).  For  the  parabolice  equation  solutions,  a  Gaussian-f ield 
initialization  was  used.  As  seen  in  Fig.  3  of  the  article,  the 
agreement  between  normal  mode  and  the  implicit  finite  difference 
with  (3.26)  methods  is  strikingly  good. 


IV.  Comments 


While  the  parabolic  equation  approximation  introduces  an 
unavoidable  phase-speed  error,  the  effects  of  this  error  can  be 
mitigated  either  by  solving  the  pseudoproblem  according  to  the 
mapping  (n,  z  )**-(n*,  z* )  developed  by  Brock  et.  al.  (1977),  or  by  the 
modified  equation  (1.29)  obtained  by  DeSant.o,  Perkins,  and  Baer 
(1978) . 

Concerning  the  split-step  PFT  algorithm,  the  depth  increment 
must  be  chosen  small  enough  that  abrupt  changes  in  the  speed  of 
sound  profile  do  not  introduce  significant  errors.  This  is  not 
too  serious  in  the  water  column,  but  the  number  of  increments  can 
become  prohibitive  in  terms  of  computational  time  if  there  is  a 
sharp  change  in  the  speed  of  sound  between  the  water  column  and 
the  bottom.  Further,  this  approach  cannot  deal  with  the  change  of 
density  between  water  and  bottom.  The  various  approaches 
attempting  to  accommodate  the  density  change  do  not  appear  to  be 
very  satisfactory,  particularly  when  the  bottom  is  an  important, 
part  of  the  propagation  problem. 

The  use  of  central-difference  methods  for  treating  the  second 
partial  of  f_  with  depth  has  opened  up  an  avenue  of  approach  which 
allows  a  much  more  realistic,  and  physically  satisfying,  treatment 
of  a  fluid  bottom.  Further  extension  along  these  lines  is  to  be 
anticipated. 

The  parabolic  equation  method  is  based  on  monofrequency, 
continuous  wave  signals.  This  means  that  the  multipath  problems 
encountered  with  transient  signals  are  not  considered,  and  to  date 
it  has  not  been  possible  to  isolate  time-separated  contributions 
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at  the  receiver  via  the  parabolic  equation.  While  in  principle 
Fourier  synthesis  over  the  required  band  of  frequencies  is 
possible,  implementation  at  this  time  would  appear  to  be  extremely 
prohibitive. 

The  assumption  of  cylindrical  symmetry  is  important  and 
restrictive  whenever  an  irregular  bottom  is  a  significant  feature 
of  the  case  of  interest.  For  example,  if  propagation  over  a 
sea-mount  is  considered,  the  assumption  of  cylindrical  symmetry 
converts  the  mount  into  a  "ring"  whose  center  of  symmetry  is  the  z 
axis.  This  results  in  the  neglect  of  any  azimuthal  reflection  or 
scattering  of  the  incident  sound  field  from  the  sides  of  the 
mount,  and  can  seriously  affect  comparisons  between  the 
cylindrically  symmetric  parabolic  equation  predictions  and 
experimental  results. 
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